// g++ 2hpc_test2.cpp -sse2 -O2
# include <iostream>
# include <cstdlib>
# include <ctime>
# include <xmmintrin.h> 
using namespace std;
//funtion1
template<typename T>
void SeqMatrixMult1(int size, T** m1, T** m2, T** result)
{
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            result[i][j] = 0;
            for (int k = 0; k < size; k++) {
                result[i][j] += m1[i][k] * m2[k][j];
            }
        }
    }
}
//function2 use temp parameter
template<typename T>
void SeqMatrixMult2(int size, T** m1, T** m2, T** result)
{
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            T c = 0;
            for (int k = 0; k < size; k++) {
                c += m1[i][k] * m2[k][j];
            }
            result[i][j] = c;
        }
    }
}

//funtion3 roll over the matrix m2
template<typename T>
void Transpose(int size,T** m)
{
	for(int i=0;i<size;++i){
		for(int j=i+1;j<size;++j){
			std::swap(m[i][j],m[j][i]);
		}
	}
}
template<typename T>
void SeqMatrixMult3(int size, T** m1, T** m2, T** result)
{
    Transpose(size,m2);
	for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
			T c=0;
            for (int k = 0; k < size; k++) {
            	c += m1[i][k] * m2[j][k];
            }
            result[i][j] = c;
        }
    }
	Transpose(size,m2);
}
//funtion4
template<typename T>
void SeqMatrixMult4(int size, T** m1, T** m2, T** result)
{
    T c;
	for (int i = 0; i < size; i++) {
        for (int k = 0; k < size; k++) {       
			T c = m1[i][k]; 
            for (int j = 0; j < size; j++) {
                result[i][j] += c * m2[k][j];
            }
//            result[i][j] = c;
        }
    }
}
//template<typename T>
//void SeqMatrixMult4(int size, T** m1, T** m2, T** result);
// assume size % 2 == 0
// assume m1[i] and m2[i] are 16-byte aligned
// require SSE3 (haddpd)

void SeqMatrixMult5(int size, double** m1, double** m2, double** result)
{
    Transpose(size, m2);
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            __m128d c = _mm_setzero_pd();

            for (int k = 0; k < size; k += 2) {
                c = _mm_add_pd(c, _mm_mul_pd(_mm_load_pd(&m1[i][k]), _mm_load_pd(&m2[j][k])));
            }
            c = _mm_add_pd(c, c);
            _mm_store_sd(&result[i][j], c);
        }
    }
  //  Transpose(size, m2);
}

// assume size % 4 == 0
// assume m1[i] and m2[i] are 16-byte aligned
// require SSE3 (haddps)
/*
void SeqMatrixMult4(int size, float** m1, float** m2, float** result)
{
    Transpose(size, m2);
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            __m128 c = _mm_setzero_ps();

            for (int k = 0; k < size; k += 4) {
                c = _mm_add_ps(c, _mm_mul_ps(_mm_load_ps(&m1[i][k]), _mm_load_ps(&m2[j][k])));
            }
            c = _mm_add_ps(c, c);
            c = _mm_add_ps(c, c);
            _mm_store_ss(&result[i][j], c);
        }
    }
    Transpose(size, m2);
}
*/
//show result
template<typename T>
void show_matrix(int size,T ** M){
	for(int i=0;i<4;i++){
		for(int j = 0;j<4;j++){
			cout<<M[i][j]<<" ";
		}
		cout<<"\n";
	}
}
int main(){
	int n = 500;
	
	double **a = new double*[n];
	double **b = new double*[n];
	double **result1 = new double*[n];
	double **result2 = new double*[n];
	double **result3 = new double*[n];
	double **result4 = new double*[n];
	double **result5 = new double*[n];
	for(int i=0;i<n;i++){
		a [i]=new double[n];
		b [i]=new double[n];
		result1[i]=new double[n];
		result2[i]=new double[n];
		result3[i]=new double[n];
		result4[i]=new double[n];
		result5[i]=new double[n];
	}
	
	for(int i=0;i<n;++i)
		for(int j=0;j<n;++j){
			a[i][j]=rand()%100/10.0;
			b[i][j]=rand()%100/10.0;
			result1[i][j]=0;
			result2[i][j]=0;
			result3[i][j]=0;
			result4[i][j]=0;
			result5[i][j]=0;
//			cout <<a[i][j]<<" "<<b[i][j]<<" ";
		}

	show_matrix(n,a);
	show_matrix(n,b);
//	show_matrix(n,result);
 
	clock_t start,end;

// method 1
	start=clock();
	SeqMatrixMult1(n,a,b,result1);
	end = clock();
	cout<<end-start<<"\n";
//	method 2
	start=clock();
	SeqMatrixMult2(n,a,b,result2);
	end = clock();
	cout<<end-start<<"\n";
//	show_matrix(n,result);
//	method 3
	start=clock();
	SeqMatrixMult3(n,a,b,result3);
	end = clock();
	cout<<end-start<<"\n";

//	method 4
	start=clock();
	SeqMatrixMult4(n,a,b,result4);
	end = clock();
	cout<<end-start<<"\n";

//	show_matrix(n,a);
//	show_matrix(n,b);	
	show_matrix(n,result1);
	show_matrix(n,result2);
	show_matrix(n,result3);
	show_matrix(n,result4);
	show_matrix(n,result5);
//delete;	
	for(int i=0;i<n;++i){
		delete[]a[i];
		delete[]b[i];
		delete[]result1[i];
		delete[]result2[i];
		delete[]result3[i];
		delete[]result4[i];
		delete[]result5[i];
	}
	delete[]a;
	delete[]b;
	delete[]result1;
	delete[]result2;
	delete[]result3;
	delete[]result4;
	delete[]result5;
	
	return 0;
}
